Statistical Eigenmode Transmission over 
Jointly-Correlated MIMO Channels 

Xiqi Gao*, Bin Jiang*, Xiao Li*, Alex B. Gershman^ and Matthew R. McKay* 



Abstract 

We investigate MIMO eigenmode transmission using statistical channel state information at the transmitter. We 
consider a general jointly-correlated MIMO channel model, which does not require separable spatial correlations at 
the transmitter and receiver. For this model, we first derive a closed-form tight upper bound for the ergodic capacity, 
which reveals a simple and interesting relationship in terms of the matrix permanent of the eigenmode channel 
coupling matrix and embraces many existing results in the literature as special cases. Based on this closed-form 
and tractable upper bound expression, we then employ convex optimization techniques to develop low-complexity 
power allocation solutions involving only the channel statistics. Necessary and sufficient optimality conditions 
are derived, from which we develop an iterative water-filling algorithm with guaranteed convergence. Simulations 
demonstrate the tightness of the capacity upper bound and the near-optimal performance of the proposed low- 
complexity transmitter optimization approach. 
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Statistical eigenmode transmission, ergodic capacity, capacity bound, MIMO channel, permanents, power 
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National Mobile Communications Research Laboratory, Southeast University, Nanjing, China 
Department of Communication Systems, Technische Universitat Darmstadt, Germany 

Electronic and Computer Engineering Department, Hong Kong University of Science and Technology, Hong Kong 



2 

I. Introduction 

Multiple-input multiple-output (MIMO) wireless systems, equipped with multiple antennas at both the 
transmitter and the receiver, have attracted tremendous interest in recent years as a means of enabling 
substantially increased link capacity and reliability compared with conventional systems [1-6]. The per- 
formance of practical MIMO systems is characterized by various system parameters, such as the average 
transmit power and the transmit-receive antenna configurations, as well as various channel phenomena 
such as spatial correlation, line-of-sight components, thermal noise, interference, and Doppler effects due 
to mobility. Each of these factors has an impact on the MIMO channel capacity. 

In realistic environments, where the channel characteristics may vary significantly over time, substantial 
MIMO capacity benefits can be obtained by tracking the states of the fading channels, and using this 
channel state information (CSI) to optimally adapt the MIMO transceiver parameters. However, such 
closed-loop MIMO strategies require both the transmitter and receiver to acquire some form of CSI. Whilst 
it is reasonable to assume that the instantaneous CSI can be obtained accurately at the receiver through 
channel estimation, whether or not this information can be obtained at the transmitter depends highly 
on the application scenario. For example, for fixed or low mobility applications, the channel conditions 
vary relatively slowly, in which case instantaneous CSI can be fed to the transmitter via well-designed 
feedback channels in frequency division duplex (FDD) systems, or using the reciprocity of uplink and 
downlink in time division duplex (TDD) systems. However, as the mobility and hence the fading rate 
increases, obtaining accurate instantaneous CSI at the transmitter becomes much more difficult. For such a 
scenario, it is reasonable to exploit statistical CSI at the transmitter. The motivation for this approach stems 
from the fact that the channel statistics vary over much larger time scales than the instantaneous channel 
gains, and the uplink and downlink statistics are usually reciprocal in both FDD and TDD systems [7, 8]. 
Therefore, the statistical information can be easily obtained by exploiting reciprocity, or by employing 
feedback channels with significantly lower bandwidth compared with instantaneous CSI feedback systems. 
In addition, transceiver designs based on statistical information are typically more robust to imperfections, 
such as delays, in the feedback channel. 

Capacity analysis and transceiver designs using the statistical CSI at the transmitter are highly dependant 
on the channel modeling. The most common approach has been to adopt the popular Kronecker model 
[9-19], where the correlation between the fading of two distinct antenna pairs is the product of the 
corresponding transmit and receive correlations [9,20]. The primary advantage of this separable model 
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is that it is analytically friendly, however various measurement campaigns have demonstrated that it can 
have deficiencies in practice [21,22]. To overcome these deficiencies, more generalized channel models 
have been proposed, including the virtual channel representation [23,24], the unitary-independent-unitary 
(UIU) model [25,26], and Weichselberger's model [21]. In contrast to the Kronecker model, these are 
jointly-correlated channel models which not only account for the correlation at both link ends, but also 
characterize their mutual dependence. 

Under various assumptions on the system configuration and channel model, several important works 
have been reported on the transmitter optimization problem using the statistical CSI at the transmitter 
in recent years. In particular, for multi-input and single-output (MISO) wireless channels with correlated 
Rayleigh or uncorrelated Rician entries, it was shown in [27] that the capacity-achieving strategy is 
to send independent data streams in the directions defined by the dominant eigenvectors of the transmit 
correlation matrix. This result was extended to Rayleigh fading MIMO channels with Kronecker correlation 
structure in [15-18], to uncorrelated Rician MEMO channels in [28,29], to the UIU model in [25,26], 
and to the virtual channel representation in [24]. These prior contributions have also considered the 
task of optimally allocating power across the transmit eigendirections (i.e. defining the eigenvalues of the 
optimal transmit covariance matrix), for maximizing capacity. However, in most cases, the power allocation 
problem has been tackled by optimizing the exact ergodic capacity expression, and this approach has led 
to computationally-involved numerical optimization procedures. For example, see [30,31] for Kronecker 
channels, and [26] for jointly-correlated UIU channels. In these contributions, iterative power allocation 
approaches were presented which involved numerical averaging over channel samples for each iteration 
of the algorithm. 

In this paper, we investigate the statistical eigenmode transmission (SET) over a general jointly- 
correlated MIMO channel, using the statistical CSI at the transmitter. Our idea is to first derive a closed- 
form tight upper bound for the ergodic capacity of the general jointly-correlated MIMO channel model. 
This upper bound expression reveals a simple and interesting relationship in terms of matrix permanents, 
and embraces many existing results in the literature as special cases, such as those presented for Kronecker 
channels in [10, 13, 14, 17, 18]. Based on this closed-form and tractable upper bound expression, we then 
employ convex optimization techniques to develop low-complexity power allocation solutions in terms 
of only the channel statistics. We derive necessary and sufficient optimality conditions, and propose a 
simple computation algorithm, inspired by the iterative water-filling techniques presented previously for 
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transmitter optimization of multiuser systems [32,33], which is shown to converge within only a few 
iterations. Numerical simulations demonstrate the tightness of the capacity upper bound and the near- 
optimal performance of the proposed low-complexity transmitter optimization approach, i.e., suffering 
negligible loss with respect to the ergodic capacity of the jointly-correlated MIMO channel. 

A. Notation 

The following notation is adopted throughout the paper: Upper (lower) bold-face letters are used to 
denote matrices (column vectors); in some cases, where it is not clear, we will employ subscripts to 
emphasize dimensionality. The N x N identity matrix is denoted by Ijv, the all-zero matrix is denoted by 
0, and the all-one matrix is denoted by 1. The superscripts (-) H , (■) T , and (■)* stand for the conjugate- 
transpose, transpose, and conjugate operations, respectively. We employ E{-} to denote expectation with 
respect to all random variables within the brackets, and use A B to denote the Hadamard product 
of the two matrices A and B. We use [A] k i or the lower-case representation a k: i to denote the (k,l)-th 
entry of the matrix A, and a k denotes the A;-th entry of the column vector a. The operators tr(-), det(-), 
and Per(-) represent the matrix trace, determinant, and permanent, respectively, and diag(x) denotes a 
diagonal matrix with x along its main diagonal. 

We will use S% to denote the set of all size-/c permutations of the numbers {1,2,..., N}, where k < N. 
By using the notation a k £ S%, we mean that a k = («i, a 2 , . . . , a k ), oti £ {1,2,..., N} for 1 < i < k, 
and di ^ otj for 1 < i,j < k and i ^ j. We will use to denote the set of all ordered length- A; 
subsets of the numbers {1,2,..., N}. By the notation a k £ we mean that a k = (a>i, (X2, ■ ■ ■ , ak), 
oii £ {1,2,..., N} for 1 < % < k, and ot\ < «2 < • • • < ctfc- The cardinalities of the sets S% and are 

w=kj. and wBw. res P ectivel y- 

With a k and (3 k defined as above, we will use A?* to denote the sub-matrix of an M x N matrix 
A obtained by selecting the rows and columns indexed by a k and (3 k respectively. A QJV will denote the 
sub-matrix of A obtained by selecting the rows indexed by a N when M > N, and As the sub-matrix 
of A obtained by selecting the columns indexed by (3m when M < N. Also, we will use a k and p k to 
denote the sequences complementary to a k and (3 k in {1,2,..., M} and {1,2,..., iV}, respectively. As 
such, A J will represent the sub-matrix of A obtained by deleting the rows and columns indexed by a k 
and j3 k , respectively. Finally, for notational convenience, we will use A»j to represent the sub-matrix of 
A obtained by deleting its i-th row and j-th column. 
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II. Channel Model and Statistical Eigenmode Transmission 
A. Channel Model 

We consider a single-user MIMO link with N t transmit and N r receive antennas, operating over a 
frequency-flat fading channel. The ^.-dimensional complex baseband received signal vector for a single 
symbol interval can be written as 

y = Hx + n, (1) 

where x is the N t x 1 transmitted signal vector, H is the N r x N t channel matrix with (i,j)-th element 
representing the complex fading coefficient between the j-th transmit and i-th receive antenna, and n is 
the iV r x 1 zero-mean additive complex Gaussian noise vector with E {nn fi } = cr^I^. It is assumed that 
x and H satisfy the following power constraints 

E{tr (xx H )} = P, (2) 
E{tr (HH ff )} = N t N r . (3) 

We define the transmit signal to noise ratio (SNR) as p = Pjo\. If the total transmitted power P is 
equally distributed across all transmit antennas, so that E |xx H } = (P/N t ) I Nt , then p also corresponds 
to the average SNR per receive antenna. 

For the jointly-correlated MEMO channel which we consider throughout this paper, the channel matrix 
H is given by 

H = U r HUf 

= U r (D + M0H« d )Uf, (4) 

where H = D + M Hud, U 4 and U r are N t x Nt and iV r x N r deterministic unitary matrices, D is 
an N r x N t deterministic matrix with at most one nonzero element in each row and each column, M is 
an N r x N t deterministic matrix with nonnegative elements, and H iid is an N r x N t random matrix with 
elements having zero mean and independent identical distributions (i.i.d.). Note that we do not constrain 
the elements of Yl iid to be Gaussian. Without any loss of generality, we can assume that the nonzero 
elements of D are real, with indices (i, i) for 1 < i < min(A r t , N r ). Let us define 

ft = E{H0H*}. (5) 
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It is easy to show that 

f2 = DoD + M©M, (6) 
and the power constraint (3) can be rewritten as 

i=l 3=1 

From (4), the transmit and receive correlation matrices can be expressed as 

R t = E{H H H} = U t A t Uf , (8) 

R r = E{HH^} = U r A r Uf , (9) 

where A t and A r are diagonal matrices with [A t ] u = J2f=i ffilji anc ^ i^-Au = Y2f=i l^lij- This implies 
that in the channel model defined in (4), U f and U r are the eigenvector matrices of the transmit and receive 
correlation matrices, respectively. These matrices are characterized by the transmit and receive antenna 
configurations. For example, when uniform linear arrays (ULA) are employed at both the transmitter and 
receiver, it is shown in [23] that the eigenvector matrices can be set to discrete Fourier transform (DFT) 
matrices. 

The statistics of H = UfHU ( characterize realistic propagation environments. From (4) and (6), we 
have 

D = E{H}, (10) 



[M] y = VHtHlii} 



["k-[D£, (ID 

where var{-} denotes variance. The matrices D and M reflect the line-of-sight (LOS) and scattering 
components of the channel, respectively. The (i,j)-th element of f2, i.e. corresponds to the average 

power of [H]y and captures the average coupling between the i-th receive eigenmode and j-th transmit 
eigenmode. For this reason, we refer to Ct as the eigenmode channel coupling matrix. It can be seen that 
the eigenvalues of the transmit and receive correlation matrices are summations of the elements of the 
matrix Cl in each column and each row, respectively. These eigenvalues are non-separable, which reflects 
the joint correlation feature of the channel. 
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The channel model described by (4) provides a general formula which embraces many existing channel 
models [9,20-26]. For example, if D = 0, M is a rank-one matrix, and Hud has Rayleigh-faded 
elements, then (4) reduces to the popular separable-correlation Kronecker model [9,20]. By allowing 
M to have arbitrary rank and fixing U t and U r to be DFT matrices, one can achieve the virtual channel 
representation for ULAs [23]. If we further allow Ut and U r to be arbitrary unitary matrices, one can obtain 
Weichselberger's channel model [21]. Moreover, by setting D = we arrive at the unitary-independent- 
unitary (UIU) model introduced in [25,26]. Our model is also related to the model in [24], where one 
LOS component was included in the virtual channel representation for the ULA MIMO channels. Here, 
we allow multiple LOS components in eigenmode to cover more general transmission links, such as those 
in distributed radio networks [34]. 

B. Statistical Eigenmode Transmission 

Throughout the paper, we assume that the receiver knows the channel perfectly, whilst the transmitter 
only has access to the statistical parameters, including U t , U r , D and M (and thus f2). Under these 
assumptions, the ergodic capacity of the MIMO channel is achieved by selecting the transmitted signal 
vector x to have zero mean and to follow a proper Gaussian distribution [1]. Let the covariance matrix 
of x be E {xx 1 } = (P/N t ) Q. Then the power constraint on x can be rewritten as tr (Q) = N t , and the 
ergodic capacity is given by 

C= max E{logdet(L Vr + 7 HQH /r )}, (12) 
where Substituting (4) into (12) yields 

C= max Ejlogdet (l Nr + 7 HQH H ) ) , (13) 

QbO,tr(Q)=7V t L \ / J 

where Q = U^QU f . Let Q = UAU^, where U is the eigenvector matrix, and A = diag(Ai, A2, ■ ■ ■ , AjvJ 
is a diagonal matrix of the corresponding eigenvalues. When H has independent and symmetrically 
distributed elements, it has been shown in [26] and [24] that the optimal eigenvector matrix U for 
achieving the capacity is U = Ut, and thus Q is diagonal. In [24], it has been pointed out that this 
solution also applies when one element of H contains a LOS component. We note however, that the 
channel model given by (4) allows for multiple possible LOS components. For this more general case, 
one can arrive at the following result: 



8 

Theorem 1: The eigenvector matrix of the capacity-achieving matrix Q for the jointly-correlated channel 
(4) is given by U = U t . The ergodic capacity can therefore be expressed as 

C = max E (log det (l N + 7 H diag(A)H^) ) , (14) 

where A is an N t x 1 vector containing the eigenvalues Aj, i = 1,2,..., N t . 

The proof follows similar approaches to those used in [24, 26, 27, 29] and is therefore omitted. Theorem 1 
demonstrates that the optimal signaling directions are defined by the eigenvectors of the transmit correlation 
matrix of the MEMO channel. This agrees with and extends prior results in the literature to the more general 
channel model given by (4). For the transmitter optimization problem, the major remaining challenge is 
to determine the eigenvalues of the capacity-achieving input covariance matrix Q. This is equivalent to 
the task of optimally allocating the available transmit power budget over the optimized transmit eigen- 
directions, determined in Theorem 1. 

In general, it is very difficult to derive an exact closed-form solution for the power allocation problem. 
A major source of this difficulty is due to the complexity in evaluating tractable closed-form solutions for 
the expectation in (14). This is also the case for many other less general MIMO channel scenarios, such 
as the popular Kronecker correlation model [15, 16]. As such, the standard approach has been to develop 
numerical optimization techniques (see e.g. [30] and [31]). 

In this paper, considering the general jointly-correlated MIMO channel model, we develop a new 
approach which leads to the design of simple, robust and practical power allocation solutions. In particular, 
our approach is based on deriving a tight closed-form upper bound on the expectation in (14) which can 
then serve as an approximation to the capacity. Based on this expression, we are then able to derive 
new optimized power allocation solutions which are simple and fast to compute. These power allocation 
solutions will be shown to serve as very accurate approximations to the optimal capacity-achieving solution, 
with low computational complexity requirements. 

We note that the power-allocation problem for jointly-correlated channel scenarios has also been 
considered in [26], where necessary and sufficient conditions as well as an iterative numerical algorithm 
were proposed. One drawback of that algorithm is that for each iteration it requires numerically averaging 
certain random matrix structures involving the inverse of instantaneous realizations of the MIMO channel. 
Moreover, since the computation algorithm requires access to instantaneous MDVIO CSI, then under the 
statistical-feedback assumption, such power-allocation computations must be typically performed at the 
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receiver. In contrast, in this paper we develop more practically appealing power-allocation algorithms 
which involve only the channel statistics. As such, they are simpler and more efficient to compute, since 
they do not require random matrix averaging during the power-allocation computation. Moreover, our new 
power-allocation algorithm has the additional advantage of permitting computation at either the receiver 
or the transmitter. This extra flexibility is particularly important for various practical applications, such as 
downlink transmission where it is often desirable or necessary to restrict computations to the base station. 
We start by rewriting the ergodic capacity (14) as 

C= max (7(A), (15) 

where 

(7(A) =E jlogdet (l Nr + 7 Hdiag(A)H H )} (16) 

is the expected mutual information between the transmitted signal x and the received signal y under SET. 
Due to the concavity of the log(-) function, the mutual information (7(A) is upper bounded by 

(7(A) < C u (\) = logEjdet (l Nr + 7 Hdiag(A)H f/ )} . (17) 

Thus, the ergodic capacity is upper bounded by 

C<C U = max (7 U (A). (18) 
*>o,E2i x *= Nt 

For the case of Kronecker MEMO channels, it has been shown in [10, 13, 14, 17, 18] that such bounds are 
very tight and admit closed-form expressions by using the expansion of the determinant. 

III. Closed-Form Capacity Upper Bound using Permanents 

In this section, we derive a closed-form expression for the capacity upper bound (18) for the jointly- 
correlated MIMO channel model in (4). We also develop algorithms for its efficient computation. The 
upper bound derivation is based heavily on exploiting linear- algebraic concepts and properties of matrix 
permanents, which we introduce and develop in the sequel. 

A. Matrix Permanents: Definitions and Properties 

The permanent of a matrix is defined in a similar fashion to the determinant. The primary difference 
is that when taking the expansion over minors, all signs are positive [35-38]. The permanents of square 
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matrices have been thoroughly investigated in linear algebra and various applied fields. The permanents 
of M x N matrices with M < N have also been defined and investigated [35]. In this paper, to facilitate 
our capacity upper bound derivation we find it necessary to extend the definition of permanents to allow 
for arbitrary M and N, and provide their useful properties. 

Definition 1: For an M x N matrix A, the permanent is defined as 



Per (A) 



M 



E YlOi.at, M < N 

N 

e n«ft„ m>n, 



(19) 



where a iy j denotes the (i, j)-th element of A. 

From this definition, one can easily establish a number of important properties of the matrix permanent, 
as given in the following lemma. These properties will be useful in subsequent derivations. 

Lemma 1: Let A be an M x N matrix, a an M x 1 vector, b an N x 1 vector, and \i a scale constant. 
Then 



Per (A) = Per(A Q 

M 

Per (a) = 



i=l 



Per(diag(a)) = det(diag(a)) 
Per(^A) = /i min(M ' Ar) Per(A) 
Per(diag(a)A) = det(diag(a))Per(A), M < N 
Per(Adiag(b)) = det(diag(b))Per(A), M > N. 



(20) 

(21) 

(22) 
(23) 
(24) 
(25) 



For an M x N matrix with M < N, there exists an analogy between the matrix permanent and the 
Laplace expansion of the determinant [35,39]. The following lemma gives the straightforward extension 
of this result for arbitrary 1 M and N. 

Lemma 2: Let A be an M x N matrix. Then 



Per (A) 



E Per ( A"* ) Per ( A J ) , M < N 

M > N, 



E Per (A?) Per ( 



(26) 



'The result for the case M > N is obtained by employing (20), and following the same steps as used in the derivation for the case 
M < N, given in [35]. 
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where a k G «s£ ) and $ k G with 1 < fc < min(Af,7V). Note that for the case k = mm(M,N), 



Per J = 1 and Per I A ? ) = L 

For the special case k — 1, (26) can be re-expressed as follows 

iV 

EflijPerfAij), M<iV 
Per(A) = { j =' (27) 
£o ili Per(A iii ), M > N, 

3=1 

where 1 < i < min(M, 7V). This is analogous to the cofactor expansion of the determinant [39]. With 
k = min(M, N), (26) simplifies to 



£ Per (A* J, M<iV 



Per(A) = < 



- c( M ) 

CTMe5jv (28) 



£ Per (A*"), M > N. 



The following two key lemmas are particularly important for deriving the closed-form capacity upper 
bound in the sequel. 

Lemma 3: Let A be an M x JV matrix. Then 

Per([I M A]) = Per([I Ar A T ] 

min(JV,M) 

= E E Per ( A " fe ) 



min(A r ,M) 

= E E Per ( A 4-)> (29) 

where Per(A afc ) = 1 and Per (As ) = 1 when k — 0. 

A proof is provided in Appendix I. The values of Per([I/i/ A]) and Per([Lv A T ]) in Lemma 3 will be 
called extended permanents of A, which we denote as 

Per(A) = Per([I A/ A]) = Per([Lv A T ]) . (30) 

Lemma 4: For an N x N random matrix X with independent elements, suppose that there exists at 
most one non-zero element in each row of X = E {X}. Then we have 

E{det(X)det (X H )} = Per (3), (31) 
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where 3 = E{X©X*}. 

A proof is provided in Appendix EE. For the special case where all elements of X are independent and 
identically distributed with zero mean and unit variance, we have that E = E{X©X*} = Inxn and 
E {det (X) det (X^)} = Per (l NxN ) = N\. This agrees with prior results in [10, 17, 18,40]. 

The following conjecture is useful when dealing with the optimal power allocation problem in Section 
IV. 

Conjecture 1: Let A be an M x A r matrix of non-negative elements. Then /(A) = log Per(Adiag(A)) 
is concave on V N = {A | Per(Adiag(A)) > 0, and \ > 0, 1 < i < N}. 

For the general case with arbitrary M and N, the formal proof of this result is not available at this 
stage. In Appendix EH, we provide proofs for several special cases, which lend support to the validity of 
this conjecture. 

B. Capacity Upper Bound 

Armed with the general results of the preceding subsection, we can now derive a closed-form expression 
for the upper bound on the ergodic capacity. 

Theorem 2: The ergodic capacity in (14) is upper bounded by 

C<C U = max logCJA), (32) 

\>0,l T \=N t 

where 

C u (\) = log Per ( 7 ft diag(A)) . (33) 
Proof: We start by writing the upper bound for the expected mutual information under SET in (17) 

as 

a(A) = log£(A) (34) 

where 

E(\) = E {det (V + 7 H diag(A)H^) } . (35) 
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By using the characteristic polynomial expansion of the determinant, as well as the Cauchy-Binet formula 
for the determinant of a product matrix [39], we have 



■ mm(N t ,N r ) , „ 

E{\) = ^\ ^ E det((fidiag(A)H^)7 



k =0 A c cW 



min(Ar t ,JV r ) 



E 



E 7*E E E ^t(Hj)de t (diag(A)f ; )det((H")^ 



k=0 Ae5« 4g5^ 



mm(N t ,N r ) 



E 7* E E det (diag(A)^)K{det (Hg)det ((H^)J)} . (36) 



fc=0 a c^O' e^W 



/ ~ \ <^fc / ~ \ 0k 

Let us denote X = H ) . . Then, X H = (H H ) , and it is easily found that E {X X*} = Q^ k . The 

\ / di, \ /at P fc 



matrix X satisfies the conditions in Lemma 4. Thus, we have 



E {det(Hj)det((H-)*)} = Per( S2 |). (37) 

Substituting (37) into (36) and using the properties of the permanents in Lemma 1, as well as (28) and 
Lemma 3, we find that 

min(N t ,N r ) 



e w= E ^ fc E E det ( dia s( A )|) Per H 



k=0 ~ c o(fc) Z c Ak) 

a k^N r Pk£& Nt 

min(N t ,N r ) 



E 7 fe E E Per ( ( " diag(A)) S) 



fc=0 a e<jM ft C?< fc ' 

min(Nt,N r ) 



7 fe E Per ((" dia g( A )r 



fc=0 a c<;W 

min(jV t ,7V r ) 



]T ^Per(( 7 ndiag(A)r 
Per ( 7 Odiag(A)). (38) 



Substituting (38) into (34) and using (18) complete the proof. □ 
From the above theorem, we see that the upper bound on capacity is completely determined by the 
average SNR p (= jNt) and the eigenmode channel coupling matrix ft. This bound is particularly useful, 
since we may now apply (33) to maximize C U (X) with the respect to A (i.e. address the power allocation 
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problem), without the need for performing Monte-Carlo averaging over random realizations of the MIMO 
channel matrix. 

It is interesting to consider the special case when D = and M = ab T , where a and b are N r x 1 and 
N t X 1 real vectors. In this case, the jointly-correlated channel model considered in this paper reduces to 
the popular Kronecker correlation model. Defining A, = a a and X t = b b, the eigenmode channel 
coupling matrix can then be expressed as ft = X r Xj = diag (A r ) lN r xN t diag (X t ), and (33) reduces to 

a(A) = logPer( 7 Qdiag(A)) 

min(N t ,N r ) 

= log 7* E E det (diag(A)|)Per ((diag (A r ) l Mt diag (A t ))£) 

min(N t ,N r ) 

= log Yl E E det(diag(A r )|)det(diag(A0A t )|). (39) 

Equation (39) is equivalent to the upper bounds presented previously for Kronecker-correlated channels 
in [17, 18]. Moreover, for the special case A = 1 (i.e. the case of equal-power allocation), (39) reduces 
further to the capacity upper bound presented in [10]. 

C. Efficient Computation Algorithms 

To evaluate the closed-form capacity upper bound expression given by (32) and (33), we must evaluate 
the extended permanent of the matrix 7f2diag(A). Clearly, when the size of the matrix is small, this can 
be done by simply expressing the extended permanent as a conventional permanent via (30), and then 
either directly applying Definition 1, or using the Laplace expansion in Lemma 2. However, in both cases, 
as the size of the matrix grows, the computational complexity increases significantly, and more efficient 
methods are needed. To see this, consider the task of evaluating the permanent of a general M x N 
matrix A. The complexity associated with computing matrix permanents is conventionally measured in 
terms of the number of the required multiplications. Adopting this measure, the number of multiplications 
required for evaluating the matrix permanent using Definition 1 and the Laplace expansion (e.g. via (27)) 
are t^z^t an d Y^k=i (■a-iey. ' res P ec ti ve ly> where m = min(M, N) and n = max(M, N). Clearly, as the 
matrix dimensions increase, the computational complexity increases exponentially. For this reason, it is 
necessary to investigate more efficient computational algorithms. 

The best-known algorithm for computing the matrix permanent of arbitrary dimensions is due to Ryser 
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[36] 2 , who showed that the permanent of the M x N matrix A (with M < N) can be evaluated via 

M M 



Per(A) = £(-l) M - fe ^_t E IN A ^ W 

fc=0 A c o(*0 1=1 

where C) = and rj(-) represents the sum of the elements in the z-th row of the matrix argument. 

j i. i). 

A similar formula also exists for the case M > N. This algorithm requires m + (m — 1) Yl™=i 
multiplications, with m and n defined as above. 

In our case, we are interested in computing the extended permanent Per (VI) in (33), i.e. the permanent 
of [Ijv r VI] or [Ijv t Vl T ], where O = 70diag(A). By directly computing this quantity based on Definition 1, 
the Laplace expansion, or Ryser's method, the number of required multiplications is ^ jVmin ~ 1 ^/ Vmi ', 1+Afmax - >! , 

''max ■ 

£S n_1 feag^fel and N ™ + " !) C k Nmin+N _, respectively, where N min = mm(N t , 

N r ) and iV max = max(N t , N r ). For practical values of N r and N t , these complexities can be quite high. 
As such, we are motivated to establish new and more efficient methods for computing the extended 
permanent, which we now consider. 

Let us define the following auxiliary matrix 

fl(z) = l NrXNt + zVt. (41) 

We will establish new efficient computation algorithms for Per(f2) based on the following key result. 
Lemma 5: Let Per(n(z)) = E£f rA) l*kZ k ''■ Then 

min(JV r ,JVt) 

Per(fi)= ^c k , (42) 

k=0 

where c k = \N t - N r \\/(m&x(N r , N t ) - k)\. 

A proof is presented in Appendix IV. This result shows that the extended permanent Per (VI) can be 
calculated directly from the polynomial expansion of Per(0(z)). Considering the case N r < N t , from 
Definition 1 in (19), Laplace expansion (27) and Ryser's expression (40), we have the following three 
formulas for Per (Vl(z)): 

N r 

Per(to{z))= E H(l+u iiai z), (43) 



qNr i = l 

5 JVt 



2 For the case of square matrices, further improvements have also been proposed [37]. 
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Nt 

Per(n(z)) = + w^Per (44) 



N r N r 



PeT(Cl(z)) = J2(-^ k C^ In&WJ. (45) 



It is convenient to re-express (45) by letting 



j(i+ ri ({i ak )z) = j2 a ^ zi > ^ 



with a ,a k = 1, such that 



i=l i=0 



l[r l (n(z) ak ) = J2k N ^a l , Qk z\ (47) 

i=l i=0 

This yields 

Nr N r 

Importantly, we find that each of the equivalent expressions (43), (44) and (48) admit simple and efficient 
recursive algorithms for calculating the coefficients of z. To demonstrate this, consider (43). Let b k (z) = 
Y\i =1 (1 + Cji. ai z) = 1 + Y.n=i bk, n z n , where k = 1, 2, • • • N r . Then, b k +i(z) = b k (z) (l + Cb k+1 ^ k+t z) for 
1 < k < N r — 1, and therefore the coefficients of z can be evaluated recursively via 



&k+l,a k+1 + b k ,i, n = 1 

Wi,n = ^ W)fc + i, au , 6fc, n -i + h,n, 2 <n < k (49) 

,6fc,fc, n = & + l. 



n — \ UJ k+ i : a k+1 b k!n - 



This result, combined with Lemma 5, presents an efficient algorithm for computing the extended permanent 
Per (CI). In a similar manner, efficient computational algorithms can also be easily obtained based on (44) 
and (48). We omit the specific details of these. For arbitrary N t and iV r , with N m i n and N max defined as 
above, the number of required multiplications for the three polynomial-based computation algorithms are 
^fr'JV' E^T" 1 ( ^T fc)J l7 x! and tf* + £ff C N , respectively. 

Fig. 1 presents the number of required multiplications for evaluating C U (X) based on the three polynomial- 
based computation algorithms, for various antenna configurations of the form N = N t = N r . The number 
of required multiplications for calculating C U (X) by directly using Definition 1, Laplace expansion and 
Ryser's formula are also shown for comparison. We clearly see that the polynomial-based algorithms have 



17 



significantly reduced computational complexity compared with the direct methods; in many cases yielding 
orders of magnitude improvements. Of the polynomial-based algorithms, the Laplace expansion gives the 
least complexity for N < 5, whereas the Ryser-based formula is most efficient for N > 5. 

IV. Optimal Power Allocation with the Capacity Bound 
A. Asymptotic Optimality at Low and High SNR 

Based on the tight closed-form capacity upper bound in Theorem 2, we can now address the transmitter 
power allocation optimization problem by dealing with only the eigenmode channel coupling matrix Q 
and the transmit SNR p (= ^N t ). The optimal solution for maximizing the upper bound will then serve 
as an approximation to the optimal capacity-achieving power allocation solution. Our numerical results 
will confirm the accuracy of this approximation. 

The power allocation optimization problem can be formulated as follows 

max C u (\) (50) 
subject to A>0, l T \ = N t . (51) 

Before dealing with this problem in its most generality, we briefly check the asymptotic optimality of our 
approach at low and high SNR. For arbitrary SNRs, we will then develop optimality conditions and an 
iterative numerical algorithm in the framework of convex optimization. 
For low SNRs, C U {X) can be expressed as 

/ N t 

C\(A)=log l+7^r l A. t + 0( 7 5 

N t 

= 7$>A i + 0( 7 2 ), ^ 



i=l 



where Tj = X^it^b*- Without any loss of generality, assume that T\ = r 2 = . . . = r t > r z+1 > . . . > r Nt . 
Maximizing the first-order (in 7) term in (52) subject to the constraint (51) gives the following power- 
allocation policy 

tN^tei-l,...,, (53) 
[0, fori = Z + l, ...,N t . 

This means that beamforming along the strongest transmit eigenmodes (specified by the channel coupling 
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matrix Q) is optimal in the low SNR case. 



For high SNRs, with N t < N r , we have 



C u (A)-^log(Per( 7 fidiag(A))) 



logPer( 7 f2) + logdet (diag(A)) 



(54) 



which is maximized by the following power allocation policy 



Ai = l, i = l,...,JV t 



(55) 



i.e. equal-power allocation over the transmit eigenmodes. These low and high SNR power allocation 
policies, derived based on the capacity upper bound, coincide exactly with the optimal capacity-achieving 
power allocation policies for the low and high SNR regimes, considered previously in [24,26]. 

B. Optimality Conditions for Arbitrary SNRs 

We now address the general case with arbitrary SNRs. To this end, let Ai > 0, A 2 > 0, and < 9 < 1. 
Then, using Conjecture 1, we can write 



where A x = [li x jv r Af] T and A 2 = [li X N r Ag] 1 '. Therefore, the function C U (X) is concave on the space 
of nonnegative A, and the optimization problem given by (50) and (51) is a concave optimization problem. 
As such, there exists only one local optimal solution, which is also a global solution. This solution could 
be evaluated by employing standard convex optimization algorithms, such as interior point methods [41]. 

Since the problem is concave, we can derive necessary and sufficient conditions for the optimal solution 
using the Karush-Kuhn-Tucker (KKT) conditions. To this end, let /i. = /x 2 , ^N t ] T an d v be the 
Lagrange multipliers for the inequality constraint A > and the equality constraint 1 T A = N t respectively. 
Then the KKT conditions satisfied by the optimal A can be expressed as 



C U (9X 1 + (1 - 0)A 2 ) = logPer^[I Wr7 Odiag(A 1 )] + (1 - 9)[I Nr 7 Odiag(A 2 )]) 



= logPer^I^T^diagCAx) + (1 - 9)[I Nr 7 fi]diag(A 2 )) 

> 9 log Per ([I Nr T fi]diag(Ai)) + (1 - 9) log Per ([l Nr 7 0]diag(A 2 ) 

= 9C U (\ 1 ) + (1-9)C U (\ 2 ), 



) 



(56) 



dC u (X) 



+ fii + v = 



(57) 
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A>0, l T X = N t , n>0, ^Ai = 0, (58) 

where 9C g£ denotes the partial derivative of C U (X) with respect to Aj, for 1 < i < N t . From (33), these 
derivatives can be written as 

dCuW 1 dE(X) 
dXi E(X) d\ ' K ' 

where E(X) = Per ( 7 f2 diag(A)). To evaluate the remaining derivatives in (59) it is useful to apply the 
Laplace expansion property of permanents, given by (27), to express E(X) as follows 

E{X)=p{\ (i) ) + \ i q{\ { i ) ), (60) 

where 

p(A(i)) = Per (70 ( i)diag(A ( i))) , (61) 

Nt 

?( A ») = X^T^.fe ( 7 n$diag(A (i) ) 

i=i 

= Per ( 7 n diag(A,)) - Per ( 7 n (i) diag(A (i) )) , (62) 

j denotes the (j, z)-th element of 17, f2^) denotes the sub-matrix of fi obtained by deleting the i-th 
column, fi^-j denotes the sub-matrix of f2 obtained by deleting the j-th row and i-th column, A(j) denotes 
the (Nt — 1) x 1 vector obtained by deleting the z-th element of A, and Aj denotes the Nt x 1 vector 
obtained by replacing the i-th element of A by unity. Therefore, (59) becomes 



dC u (X) g(A (i) ) 



d\ p(A(i)) + \q(\i)) ' 
Substituting (63) into (57) and eliminating the slack variable /x, the KKT conditions become 



(63) 



1 T A = N t , (65) 

where (a) + = max{0, a} and v = l/u. 

In summary, we have the following theorem. 

Theorem 3: The expected mutual information upper bound C U (X) is concave with respect to A, and 
the necessary and sufficient conditions for optimal power allocation are given by (64), where v is chosen 
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to satisfy the power constraint in (65). 

Note that when the eigenmode channel coupling matrix fl is square and diagonal 3 , we have 



and the conditions in (64) simplify to 



g(A(j)) = w iti p(\fl) , (66) 



\— [v — ] ■ (67) 



/.I 



This is the same formula as the water- filling solution when the transmitter has instantaneous CSI [1], and 
one can easily obtain the optimal power allocation via the water-filling algorithm. However, in the general 
case of an arbitrary eigenmode channel coupling matrix, the solution can not be obtained as easily and 
numerical approaches are required. 



C. Iterative Water-Filling Algorithm 

In this section, we propose a simple iterative water-filling algorithm (IWFA) for evaluating the optimal 
power allocation policy which satisfies (64). Our algorithm is based on observing that the right-hand side 
of (64) is independent of Aj, and is motivated by the IWFA methods proposed in [32, 33] for transmitter 
optimization of multiuser systems with instantaneous CSI known to the transmitters. Simulation results, 
given in Section V, show that this approach works very well and is highly efficient; typically converging 
after only a few iterations, with the first iteration achieving near-optimal performance. The proposed 
algorithm includes the following steps: 

(1) Initialize A = 1, C U (X°) = logPer( 7 ri), and k = 0. 

(2) Calculate p(A^) = Per( 7 ft (i) diag(Af )) and q(X k (i) ) = Per ( 7 diag(A 4 fe )) - Per( 7 ft (i) 

diag(A^)),z = l,2,...,iV t . 

p(\ k ) 

(3) Calculate X k+1 — (y — ttfM + ' * — 1>2, N t , via the conventional water-filling algorithm with 
power constraint Ylf=i = Nt. 

(4) Calculate C it (A fe+1 ) = log Per( 7 fi diag(A fc+1 )) . 

(5) If C u (X k+l ) < C u (\ k ), set A fc+1 := ^A fc+1 + ^\ k , and recalculate C u (\ k+l ). 

(6) Set k := k + 1 and return to Step 2 until the algorithm converges or the iteration number is 
equal to a predefined value. 

3 In this special case, the MIMO channel is essentially reduced to a set of non-interfering scalar subchannels. 
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Here, A fc stands for the value of A in the k-th iteration. In Step 1 in the first iteration, A is initialized to 1, 
i.e., to the equal power allocation. Note, however, that A could also be initialized in a different way. For 
example, in practice it is reasonable to suppose that the channel statistics change smoothly from frame 
to frame, where a more appropriate starting point for any given frame would be the optimal value of A 
from the previous frame. This could speed up the convergence of the IWFA. In Step 3, the conventional 
water-filling algorithm is performed with the required variables p(\u)) and q{\{\) calculated in Step 2. 
Following the calculation of the C U (X) in Step 4, Step 5 is performed to guarantee the convergence of the 
iterative procedure. We discuss this issue in detail below. In Step 6, the convergence of the algorithm can 
be determined by checking whether \C u (X k+1 ) — C u (X k )\ (or ||A fe+1 — A fe ||) is less than some predefined 
tolerance. 

Theorem 4: The IWFA for optimal power allocation converges to the capacity upper bound C u . 
Proof: In order to verify the convergence of our proposed IWFA for optimal power allocation, we 
define the following function for a given X k : 

1 Nt 

= w^ log (p(A « } + A ^ (A » } ) • (68) 

* i=l 

It can be seen that C U (X) is a concave function with respect to A. The water- filling solution in Step 3 
of the IWFA is exactly equal to the solution of maximizing C U (X), for a given X k subject to the power 
constraint 1 T A = N t . Therefore, with the A fc+1 resulting from Step 3 of the IWFA, we have 

C u (X k+1 ) > C u (X k ) = C u (X k ). (69) 

From the concavity of C U (X), it can be shown that the following relation holds: 

a(A fc+1 ) < C u (i-A fc+1 + . (70) 

Combining (69) and (70) yields 

C u (X k ) < C u ^X k+1 + . (71) 

Therefore, after Step 5 of the IWFA, we have that C u (X k+l ) > C u (X k ). This, along with the fact that the 
problem (50)-(51) is convex, completes the proof. □ 
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Notice that the relation (71) suggests, mathematically, to update A with (-^-A fc+1 + ^j^-X k ) in the k-th 
iteration of the IWFA, whereas the KKT conditions (64) suggest a more intuitive interpretation based 
on the water-filling principle. In our proposed IWFA, we update A with the water-filling solution if the 
resulting C u (\) is increased. This allows very fast convergence, as we demonstrate through simulations in 
the following section. To guarantee the convergence, we use (^X k+1 + ^j^-X k ) to replace the water-filling 
solution when the resulting C U (X) is not increasing in each iteration. 

V. Simulations 

In this section, we present numerical results to evaluate the tightness of the capacity bound, and to 
demonstrate the efficiency and performance of the proposed transmitter optimization approach under SET. 
We consider a MIMO system with five transmit and five receive antennas, and present results for both the 
jointly-correlated MIMO channel model and the Kronecker-correlation model. For the jointly-correlated 
channel, we adopt the same channel parameters as used in [24], where D = and Q has the following 
structure 
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For the Kronecker channel, we adopt the constant-correlation model for constructing the transmit and 
receive correlation matrices [10]. An N x N constant-correlation matrix is given by 



®n (a) 



al 



NxN 



+ (1 - a) I 



N; 



(73) 



where a E [0, 1] is the correlation coefficient. We set the transmit and receive correlation coefficients to 
be a t = 0.4 and a r = 0.6 respectively. 

Fig. 2 compares our closed-form ergodic mutual information upper bound (18) with Monte-Carlo 
simulated exact curves based on (16), for the case A = 1 (equal-power allocation). Results are shown for 
both the jointly-correlated channel and the Kronecker channel, with the above settings. We see that the 
upper bound is rather tight for both channel models, especially for low to moderate SNRs (eg. < 8 dB). 
Moreover, we see that the bound for the Kronecker model is slightly tighter than for the more general 
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jointly-correlated model. Interestingly, we will show that, despite this difference in tightness, the low- 
complexity power allocation policies derived based on these bounds perform near-optimally for both the 
Kronecker and jointly-correlated channel models. 

Fig. 3 and Fig. 4 present the ergodic mutual information achieved by the SET approach employing the 
proposed IWFA (derived based on our closed-form upper bound), in the jointly-correlated and Kronecker 
channel scenarios, respectively. For comparison, the exact ergodic capacity curves are also shown, which 
were obtained by numerically evaluating (14) using a constrained optimization function of the Matlab 
optimization toolbox. The ergodic mutual information achieved with equal power allocation (55) and 
beamforming (53) are also shown for further comparison. We clearly see that, for both channel models, 
the proposed SET approach performs near-optimally, suffering almost negligible loss compared with the 
true channel capacity. Furthermore, we see that equal power allocation and beamforming are optimal in 
the high and low SNR regimes, respectively, which agrees with our analytical conclusions put forth in 
Section IV-A. The capacity upper bound curve is also shown on the figures, and once again is seen to be 
tight. 

Fig. 5 and Fig. 6 demonstrate the convergence of the proposed IWFA for optimal power allocation 
in the jointly-correlated and Kronecker channel scenarios, respectively. Here, the SNR p was set to 10 
dB, and in all cases the algorithm was initialized using A = 1. These figures show the evolution of 
the eigenvalues Aj, i = 1, ... ,5, and the capacity bound C U (X) for each iteration. From these results, 
we see that the proposed IWFA converges after only a few iterations, with the first iteration achieving 
near-optimal performance in all cases. 

VI. Conclusions 

We have investigated statistical eigenmode transmission over a general jointly-correlated MIMO channel. 
For this channel, we derived a tight closed-form upper bound for the ergodic capacity, which reveals a 
simple and interesting relationship in terms of matrix permanents of the eigenvalue coupling matrix, and 
embraces many existing results in the literature as special cases. Based on this expression, we proposed 
and investigated new power allocation policies in the framework of convex optimization. In particular, 
we obtained necessary and sufficient optimality conditions, and developed an efficient iterative water- 
filling algorithm with guaranteed convergence. The tightness of the capacity bound and the performance 
of our novel low-complexity transmitter optimization approach was confirmed through simulations. Our 
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approach was shown to suffer near-negligible loss compared with the ergodic capacity of the jointly- 
correlated MEMO channel. 

Appendix I 
Proof of Lemma 3 

Let I = [1 M Ojuxjv] and A = [0 Mx m A]. From the definition of the permanents, we have 

Per([I M A]) = Per(A + I) 

M 

= E H^m,0 m +im,p m ), (74) 
Pm€S*< +n m=l 

where i mn and a m n denote the (m, n)-th elements of I and A respectively. Note that the following 
identity holds: 

M M k M-k 

[ {x m + y m ) = E I ] [ y*' m i (75) 



m = l k=0 - ,- c (fc)m=l m=l 



where a 2 , a M-k) e Slj k is the sequence complementary to in {1,2,..., M}. Hence 



M 



( k M-k \ 



Per([I M A]) = ]T 



(76) 



e n n a m tP „ 

It can be seen that i_ a > p , = 6((3 a > -a'j, where £(■) is the Kronecker delta operator, and Um=i &* m ,p am ^ 
only if (3 am > M and k < min(M, N). Therefore, we have 

min(M,A r ) k 

Per([I A/ A])= E E E II 

k=0 ~ r- cCO B, cS k "1=1 
min(M,N) 

= E E E Per (aJ) , (77) 



where Per ^A^ fc j = 1 when k = 0. Using (28), we have 



min(M,7V) 

Per([I M A])= J2 E Per (A" fc ) 



k=0 i, eS W 



min(M,7V) 

E E Per ( A 4)- (78) 
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Through a similar procedure, one can obtain that 



min(M,7V) 

Per([l^A T ])= E Per (A^ 



k=0 , c «jW 



min(M,iV) 

E E Per ( A /U- ( 79 ) 



fc=0 

This completes the proof. 

Appendix II 
Proof of Lemma 4 

From the definition of the determinant, we have 

E{det(X)det(X")}= E £(-ir< 4 ">(-l)^E \fl x ^b\ > ( 80 ) 

where cr(«iv) denotes the number of inversions in the permutation oln from the normal order 1, 2, . . . , N, 
and Xij is the (i, j)-th element of X. Since the rows of X are independent, we have 

{N ~\ N 

j=l J i=l 

Since the elements in each row are independent and there is only one possible non-zero mean element in 
each row, we have 

E { x i,<*i x ipi} = - (82) 

where £jj is the (z, j)-th element of S. Substituting (82) into (81) and then into (80) yield 

TV 

E{det(X)det(X")} = ^ J] = Per (E) . (83) 

& N dS% »=1 

This completes the proof. 

Appendix III 

Proof of the concavity of /(A) = logPer(Adiag(A)) in several cases 

Case 1: M > N. In this case, we have that /(A) = logPer(A) + logdet(diag(A)). The concavity of 
/(A) comes from that of logdet(diag(A)) [41]. 
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Case 2: M = 1 and N > 1. In this case, A is a row vector, and we have that /(A) = log(AA). The 
concavity of /(A) comes from that of the log function. 

Case 3: M = 2 and N > 2. In this case, we will first show that the following inequality holds: 

g 2 (X 1 + A 2 ) > g 2 (X 1 ) + g 2 (X 2 ) 
gi(Xi + X 2 ) ~ 0i (Ai) 0i(A 2 )' 

where 01(A) = li x2 Adiag(A) and 02(A) = Per(Adiag(A)). Then we will prove the concavity of /(A) 
from (84). 

Since 01(A) and 02(A) are positive on V N , the inequality (84) holds if and only if the following 
inequality does: 

£(Ai, A 2 ) = g 2 (X 1 + A 2 )0i(A 1 )^ 1 (A 2 ) - 2 (A 1 )^ 1 (A 1 + X 2 ) gi (X 2 ) - 2 (A 2 )0i(Ai + X 2 )gi(Xi) > 0. (85) 

Let A = {&[ sl 2 ] t . Then we have that 01(A) = af A + aJX and 02(A) = aJXa^X — A T diag(a! a 2 )A. 
By substituting these expressions into 0(Ai, A 2 ), we can obtain 

g(X 1 , A 2 ) = (af X 1 a 2 v X 2 - af A 2 a^Ai) 2 + (gi(X 2 )Xi - 5fi(A 1 )A 2 ) T diag(a 1 a. 2 )(gi(X 2 )X 1 - 0i(Ai)A 2 ). 

(86) 

Therefore we achieve (85) and then (84). From (84), we have 

a,(0Ai + (l-fl)A 2 ) > , 2 (Ai) 2 (A 2 ) 

g 1 (6X 1 + (l-e)X 2 ) - 0i(A x ) +l >g x {X 2 y ^ ; 

where < 9 < 1. Taking logarithm on both sides and using the concavity of the log function yields 

/(0Ai + (1 - 9)X 2 ) - 6f(X 1 ) - (1 - 9)f(X 2 ) 

> log(g 1 (9X 1 + (1 - 0)A 2 )) +9log(g 1 (X 1 )) -(1-9) log( 9l (X 2 )) > 0. (88) 

This completes the proof of the concavity of /(A). 

Case 4: A is of rank one. Let A = ab T , where a and b are vectors of M and N elements respectively. 
In this case, we have 

/(A) = logPer(l A/x7V diag(b © A)) + logdet(diag(a)) 

= log ^ Per ((lA/xwdiag(b X)) &N ) + log det(diag(a)) 

c*m£S n 

= \ogE M (b A) + log(M!) + logdet(diag(a)), (89) 
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where the function E M (\) is the M-th elementary symmetric function defined as [42] 

M 

e m {\)= II A « < 90) 

Since Em{^) is logarithmically concave, we obtain from (89) that /(A) is concave. 

Appendix IV 
Proof of Lemma 5 

We consider the case with N r < N t . The proof for the case with N r > N t is similar. From the definition 
of the permanents, we have 

N r 

Per(h(z))= Q(1+^W), (91) 



,- c N r 7 — 1 

a Nr £S N ^ 

where cDjj represents the (i,j)-th element of fi. For each product term in the above expression, the 
following relation holds: 

N r N r k 

n(i+£wo=xx e n ( 92 ) 

i=l fc=0 4e5 « i=l 

Substituting (92) into (91) yields 



Per(ft(.)) = 5> fe £ £ 



fc=0 A ,<.W A „ (zC^Vr 1 = 1 

Nr (-/V" /l) ! ^ 

= E^ (A r/_ A r r)! E E 11^ 

= E**$^Sf £ (93) 

From Lemma 3, we have the expansion of Per(fi). By comparing the resulting expansion of Per(fi) with 
(93), we complete the proof. 
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Fig. 2. Comparison of the exact ergodic mutual information and the mutual information upper bound, 
correlated channel model and the Rronecker model, with equal-power allocation (i.e. A = 1). 



Results are shown for the jointly- 
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Fig. 3. Comparison of the ergodic capacity of the jointly-correlated MIMO channel achieved by numerically solving (14), and our proposed 
iterative water-filling algorithm under SET. The capacity upper bound and the information rates achieved by equal power allocation and 
beamforming are also shown. 




Fig. 4. Comparison of the ergodic capacity of the Kronecker MIMO channel achieved by numerically solving (14), and our proposed 
iterative water-filling algorithm under SET. The capacity upper bound and the information rates achieved by equal power allocation and 
beamforming are also shown. 
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Fig. 5. Convergence of the iterative water-filling algorithm for optimal power allocation in the jointly-correlated channel. Results are shown 
for SNR = 10 dB. 




Fig. 6. Convergence of the iterative water-filling algorithm for optimal power allocation in the Kronecker channel. Results are shown for 
SNR=10 dB. 



